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ABSTRACT 


A simple rectangular finite element was developed for two-dimensional 
analysis of laminated composite materials. The rectangular laminated com- 
posite (RLC) element eliminates the need to add elements to a model simply to 
account for the material properties of various laminae. This is particularly 
advantageous for thick laminates with many lamina. Explicit integration in 
terms of generalized displacements minimizes the algebraic effort required to 
derive the element stiffnesses and the thermal load vector. A substitute 
shape function technique is used to improve the performance of the element in 
modeling bending type deformation. Results for several example problems are 
discussed. 
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INTRODUCTION 


The development of an appropriate finite element mesh Is a key step In 
successful finite element analysis* For homogeneous materials the mesh 
refinement Is dictated by geometrical considerations. The shape of a 
structure should be faithfully modeled. Also, extra mesh refinement Is 
required In regions with strong strain gradients caused by holes, cracks, or 
boundary conditions. For laminated materials the analyst must also account 
for the different material properties of the various laminae. These Ideas are 
Illustrated by the laminated composite beam shown In figure 1. Geometrical 
considerations require very few elements except close to the point where the 
load Is applied, where strain gradients are large. But since standard finite 
elements cannot account for stacking sequence effects such elements should not 
span across lamina boundaries. Hence, because of the laminated character of 
the material, the mesh should be highly refined even where the strain 
gradients are small. 

The expense of modeling each lamina individually rapidly becomes 
Intolerable as the number of laminae Increases. Reference 1 presents an 
approximate technique to reduce costs. Laminate theory is used to obtain 
effective extensional moduli for a group of laminae. Then the group of 
laminae rather than the individual lamina are modeled using finite elements. 
This approach ignores stacking sequence effects within the lamina group. 
Therefore, the flexural and flexural-extension coupling properties of the 
lamina group cannot be faithfully modeled. Reference 2 presents a hybrid 
analysis for thick laminates. In this analysis (which Is not a finite element 
®*^3lysls) the laminate Is divided Into global and local regions. The terms 
global and local refer to the detail with which the Individual lamina is 
modeled; the local region Is modeled with much greater detail than the global 



region. Conceptually, this is similar to using a finite element model with 
smaller- or higher-order elements In one region than in another. However, the 
analysis in reference 2 does not offer the inherent flexibility of the finite 
element method for modeling complicated geometries and for performing con- 
vergence checks. The objective of this paper is to introduce a new type of 
two-dimensional (l.e., plane stress or plane strain) finite element for 
analysis of laminated composites. 

The element is a four-node, bilinear, rectangular element. An ordinary 
bilinear rectangle performs poorly in modeling bending-type deformation. The 
performance can be improved by using reduced numerical integration or substi- 
tute shape functions (ref. 3). Because of the multiple laminae within the 
element, numerical integration is not appropriate. Therefore, substitute 
shape functions are used to improve the performance. Explicit integration of 
the element stiffness matrix in terms of generalized displacements minimizes 
the algebraic effort required to account for the various laminae within a 
single element. 

After describing the theoretical aspects of the element, results from 
analyses of several simple configurations are discussed. 
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NOMENCLATURE 


A area 

^11»M2*^22»^33* 

atb,C)d,a)b,C|d, 


coefficients related to laminae material properties 
generalized degrees of freedom 




F 

F 

^n 

% 

H 

K 

K 

Kn, 


m 




‘N 

Kg 

V 

N 

T 

Tn 


t 

U 


plane stress or plane strain material stiffness coefficients- 
ij - 1,3 

force vector 

transformed force vector 

force corresponding to degree of freedom n 
subvector of force vector related to normal strains 
subvector of force vector related to shear strain 
matrix used in calculation of transformation matrix 
element stiffness matrix 
transformed element stiffness matrix 

stiffness term in K, n,m = 1, nximber of degrees of freedom 

submatrix of stiffness matrix related to normal strains 

submatrix of stiffness matrix related to shear strain 

half-length of element In x-direction 

half-length of element in y-direction 

number of plies 

transformation matrix 

transformation matrix for normal strain related terms In generalized 
stiffness matrix and force vector 

transformation matrix for shear strain related terms in generalized 
stiffness matrix and force vector 

thickness in z-direction 

strain energy 
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u displacement in x-direction 

uj^,U2,U3»u^ nodal displacements in x-directlon 

V displacement in y-direction 

Vi,V2»V3,v^ nodal displacements in y-directlon 

W potential energy of external loads 

x,y rectangular Cartesian coordinates 

y-coordinate of bottom surface of i*"^ ply* Bottom ply in element is 
ply 1 . 

A vector of generalized displacements 

Ajj generalized degree of freedom n 

& vector of nodal displacements 

H strains (e^ = 

^x’®y»^xy strains 

n total potential energy 

Definitions 

generalized displacements parameters related to translation, rotation, 

and deformation of an element but not 
associated with a node 

generalized forces force associated with generalized 

displacement 

generalized stiffness matrix stiffness of an element in terms of 

generalized displacements 
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THEORY 


The RLC element is a rectangular element with four nodes. All laminae 
in the element are assumed to be orthotropic, oriented parallel to the 
x-axis, and to extend across the element width (fig. 2). The origin of the 
x-y coordinate system is at the element centroid. 

The following sections derive the RLC element stiffness matrix and the 
equivalent nodal load vector for initial thermal strains. The derivation 
begins with the presentation of general expressions for element forces and 
stiffnesses for an arbitrary finite element. Next, the particular shape 
functions used to approximate the displacements and strains are discussed. 

Then explicit expressions for the stiffness matrix and element forces due to 
thermal strains are derived. These expressions are in terms of generalized 
displacements. The final section describes how to transform the stiffness 
matrix and forces from a system of generalized displacements to one of nodal 
displacements. 

Cartesian tensor notation is used herein to express several of the 
complicated equations in compact form. In these compact equations the 
strains and refer to e^, e^, and respectively. Also, 

some parameters refer to an entire vector or matrix when there is no subscript 
(eg. F) and to a single value when there is a subscript (eg. Fj^). 
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General Expressions 

The total potential energy, n , Is given by equation (1) (ref. 1): 


n 


u + w 


t 

2 





dA - F 

n 


A 


n 


( 1 ) 


where U Is the strain energy and W Is the potential energy of the applied 
loads. Minimization of H with respect to the generalized degrees of freedom 
(d.o.f.), Ajj, yields the generalized force associated with each d.o.f.. 
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( 2 ) 


The terms In the element stiffness matrix are calculated by differentiating 
the generalized forces with respect to the d.o.f., equation (3). 
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Since linear strain-displacement relations are used In this paper, the 

9A 3A zero* Therefore, the terms in the element stiffness matrix 

n m 

are 
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nm 
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9e, 9e. 
m n 


(4) 


Shape Functions and Strain Expressions 
The technique of substitute shape functions was used to Improve the 
performance of the RLC element In modeling bending type deformation (ref. 3). 
This technique Involves using different shape functions for terms related to 
normal strains and for those related to shear strains. 
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The shape functions used in calculating terms related to normal strains 
are given by equations (5); 

u “ a + bx + cy + dxy 

(5) 

V = a + bx + cy + dxy 


The normal strains, and 


Cy, are therefore 


3u 

^x 3x 
3v 

C a — 

y 3y 


b + dy 
c + dx 


( 6 ) 


The shape functions used in calculating terms related to the shear strain 
are given in equations (7) 


u = e + fx + gy 


V = e + fx 4- gy 


The shear strain, e^^y, is therefore 

3u . 3v „ „ . 
e = -T— + -r- = g + f 
xy 3y 3x 


(7) 


( 8 ) 


Equation (8) shows that the shear strain is constant for the RLC element. 
This constant value is defined to be ”h” in equation (9). 


= h = g + f (9) 

Element Stiffness Matrix 

Identification of the relevant d.o.f. is the first step in the calcula- 
tion of the element stiffness matrix. Equation (4) shows that K is only a 
function of those d.o.f. used to define the strains. Hence, the relevant 
d.o.f. are 
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'b 

d 


d 

h 


( 10 ) 


The 5x5 element stiffness matrix can now be calculated using equations (4), 
(6), (8), (9), and (10). The non-zero terms In the stiffness matrix are 


^11 “ ^11 

^12 “ \i 

^13 ^12 

K22 * 2 l^t 

^23 “ ^Ijjt Bj_2 


44 

55 


K21 = 


K. 


31 


K32 = 


2 3 

3 V ^22 

^33 

Ki2 

K13 

^23 


(11) 


%3 “ A22 


The remaining are zero. 
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where 


^11 " £ ^11 1+1 " 
i-1 c. 

^12 “ ^ ^12^^i+l " 
i»l 

N ^ ■ 

^22 " 21 ^22^^!+! " 

1-1 

^3 " X ^^(^1+1 “ V. <^2) 

. 1-1 

®ii “iX ^u(yi+i - yi) , 

1=1 

°12 ■ ?X ■ ’'i' 

1=1 

' “u-?i; <i(.yh-yl) . 

1=1 

where 

t = element thickness 

9 l^ = half-width of element (see fig. 2) 

iy = half-height of element (see fig. 2) 

N = number of laminae 

Cf- = plane stress or plane strain stiffness coefficients, (i,j = 1»3) 
y coordinate of bottom surface of ply number ” 1 ’*, (1 = 1 ,N) 

Note that in equations ( 11 ), the only nonzero term. in K related to the 
shear strain is K55. 
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Thermal Load Vector 


Equation (2) gives the general expression for element forces. As was 
the case for the stiffness matrix, the relevant d.o.f. for thermal loads are 
given by equation (10). The strain sj to be used in equation (2) is 

Sj = Oj AT (13) 

Because the material is orthotropic, 03 = 0. The element forces 
Corresponding to each d.o.f. can now be calculated using equations (2), (6), 
(8), (9), (10), and (13). These forces are 

’'1 “ 'E (=u<4 “"Kyi+1 - yi> 

1=1 

”2 - “x Z ’^24 

1=1 

N 

’'a = + 44 4i‘Ky„.i - y^) 

F4 = 0 
F5 = 0 

Note that the force corresponding to the shear strain related d.o.f., Fc Is 
zero. 
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Transformation of Element Stiffnesses and Forces 

The preceding sections give expressions for element forces and stiff- 
nesses for a system with generalized d*o«f« A. However, to assemble element 
stiffnesses and forces into a global system of equations requires that the 
d.o.f. be nodal displacements, not generalized displacements. 

Equations (15) give the general procedure for transforming the stiffness 
matrix and force vector from one set of d.o.f. (A) to another set (6): 

■K* = T%T 

F’ = T^-F 

where T is defined by the equation: 

A = T6 

In equations (15) , K and F are in terms of generalized displacements A 
and K' and F' are in terms of nodal displacements 6. The matrix T is 
the transformation matrix. 

The first step is to calculate the transformation matrix T. Since the 
displacements u and v are approximated by different shape functions for 
terms related to normal and shear strains, the transformation matrices for 
these two types of terms must be calculated independently. 

The transformation matrix for the terms related to normal strains is 
calculated first. Equations expressing the nodal displacements u^, 

1 terms of all the element's generalized displacements are given in 

equations (16) . 
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0 1 


-i z 0 

X y 


Z Z 0 

X y 



0 0 0 


— — 
a 

11 


b 

X y X y 



0 0 0 


c 

1 -1 -11 


d 

X y X y 



0 0 0 


a 

11 11 


b 

X y X y 



0 0 0 


c 

-1 1-11 


d 

X y X y 











( 16 ) 


Equations (16) are formed using the expressions in equations (5). Equa- 
tions (16) can be solved for A« 
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0 

1 1 
X y 

0 

*xV 

0 


"1 

b 



0 


0 

V 

0 

~^y 

0 


^1 

c 


-^x 

0 

-‘x 

0 

*x 

0 

^x 

0 


“2 

d 

1 

1 

0 

-1 

0 

1 

0 

-1 

0 



a 

4l 1 
X y 

0 

1 1 
X y 

0 

*xV 

0 

‘x'y 

0 



“3 

b 


0 

-1 

y 

0 


0 


0 

-''y 


^3 
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0 

-1 

X 

0 

-1 

X 

0 

‘x 

0 

*x 


“4 
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0 

1 

0 

-1 

0 

1 

0 

-1 


^4 

A 
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r 

-‘*xV 





't, 


(17) 
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Comparison of in equation (17) and T In equation (15) shows that 

T * H Since the generalized forces and stiffnesses Involve 

only the d.o.f. b, d, c, and d, only rows 2 , 4, 7, and 8 of are 
required for the transformation matrix. Hence, the transformation matrix for 
terms related to normal strains is 



0 









1 

0 


0 -1 0 1 
0 0 


0 -1 
^ 0 


0 



1 0-1 0 1 


0 



(18) 


Now the transformation matrix for the terms related to shear strains is 
derived. As before, the first step is to express the nodal displacements in 
terms of the generalized displacements; in this case equations (7) are used. 




(19) 
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Since the number of equations exceeds the number of unknowns, equations (19) 
cannot be solved for A In terms of 5. Equations (20) give a new set of 
equations which are equivalent to equations (19) in a least squares sense 
(ref. 3). 

h'^6 = A (20) 

Solution of equations (20) yields: 



( 21 ) 


The parameters g and f are the only relevant d.o.f. These two are com- 
bined to form h (see eq. (9)). Therefore the transformation matrix for the 

T T 

shear related terms is obtained by adding rows 3 and 5 of [H H] 



( 22 ) 
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The transformation matrices % and Tg are combined to form the total 
transformation matrix, as Illustrated In equations (23). 



(23) 



The stiffness matrix K* In equation (23) was derived assuming there 
were multiple laminae within the element. If there Is only one lamina (or all 
the lamina are Identical), the stiffness matrix Is Identical to that for an 
ordinary bilinear rectangular element with reduced Integration. 
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EXAMPLE PROBLEMS 


This section discusses results from analyses of several simple problems 
using the RLC element. All the configurations consist of a laminated canti- 
levered beam (fig. 3) with either mechanical or thermal load. The canti- 
levered beam configuration was chosen because it is a severe test for 2— D 
plane stress (or plane strain) elements. Several combinations of lamina 
properties and stacking sequence were examined. Material properties for the 
various laminae are given in Table 1. All laminae were assumed to be 
isotropic. Lamina types H, I, and J have the same Young's modulus and shear 
modulus only the thermal expansion coefficients are different. Lamina type S 
has 104 of the stiffness of the other lamina types. Figure 4. shows the five 
finite element meshes used in the convergence study. Results for mechanical 
and thermal loads are discussed separately in the following sections. 

Mechanical Load 

Three laminates were examined: (H/H/H/H), (H/S/S/H) , and (S/H/H/S). The 

first two laminates have about the same flexural stiffness. The third 
laminate is much more flexible, since the outer plies are soft. 

The loading consisted of a single point load at the end of the beam. For 
a long, thin cantilevered beam (such as that in fig. 3), the tip deflection 
calculated using strength of materials, is given as: 

Pf,^ 

^ - W- (24) 

where P *» applied load 
Z » beam length 
D » flexural stiffness 

“ tip deflection calculated using strength of materials 
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Assuming equation (24) gives the correct solution, figure 5 shows the 
error in the calculated tip deflection for. the three laminates using the five 
meshes in figure 4. The open symbols show the results obtained using the RLC 
element described in the preceding sections* The solid symbols are for a 
modified RLC element and will be discussed later in this section. The error 
reduces rapidly with increased mesh refinement. The laminate (H/H/H/H) has no 
lamination effects — since all the layers are the same. As pointed out 
earlier, the element stiffness matrix is therefore identical to that for an 
ordinary bilinear rectangular element with reduced integration. Figure 5 
shows that the error for the nonhomogeneous laminates (ie. (H/S/S/H) and 
(S/H/H/S)) is comparable to that for the (H/H/H/H) laminate. Therefore, 
additional errors due to lamination within an element appear to be small. 

Much of the error in the results is due to the assumption that Cy 
within an element does not vary in the y-dlrectlon (see eq. (6)). Because of 
the Poisson effect, the upper part of the beam (which has positive should 

have a negative Sy. The lower part of the beam (which has a negative 
should have a positive gy. But within a single element gy is constant in 
the y-direction. Therefore, if there is only one element through the thick- 
ness, gy is calculated to be zero. This results in an overly stiff 
element. The magnitude of the error for a homogeneous isotropic can be esti- 
mated by examining the constitutive equations. Assuming plane stress condi- 
tions for an isotropic material, the stresses can be expressed as 



+ V£y) 


( 25 ) 


I - V 
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For a long, thin beam, a„ should be negligible; therefore, e “ -ve and 

j y X 

“ Ee„. But If z Is constrained to be zero, then 


"x 1 2 X 

1 - V 


(26) 


This results In an effective modulus of E/(l - v^). For v ®* 0.3 this 
produces a modulus which Is 10% too large. Since the flexural stiffness Is 
linearly related to the modulus, the deflection Is Inversely and linearly 
related to the modulus (see eq. (24)). Therefore, a minimum of about 10% 
error Is expected when 1 element Is used through the thickness of the beam. 
Figure 5 agrees quite well with this prediction. Of course, with two elements 
through the thickness, the spurious stiffening Is less (see fig. 5). Spurious 
stiffening due to Poisson's effect can be eliminated by artificially setting 
Vjjy = 0. When Is artificially set to zero, the element will be referred 

to as the modified RLC element. The solid symbols In figure 5 show results 
obtained with the modified RLC element. The convergence is seen to be 
extremely rapid. Further testing of the modified RLC element Is needed to 
determine when the artificial prescription of Vjjy «= 0 may lead to problems. 


Thermal Load 

Two laminates were examined: (H/H/I/I) and (H/H/J/J). In the first 

laminate, the Initial thermal strains In the top two laminae are simply the 
negative of the initial thermal strains in the bottom two laminae. The second 
laminate is like the first except that the initial thermal strain Cy Is the 
same for all four laminae. Using elementary beam theory, the end displacement 
is independent of the initial strain Sy. Hence, both beams should have the 
same end displacement. The technique for solving this problem using strength 
of materials is outlined in reference 4 and will not be discussed here. 
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Figure 6 shows the accuracy of the RLC element in calculating end 
displacement for the (H/H/I/I) laminate. The strength of materials solution 

Is assumed to be correct. The open and solid symbols are results for the 
standard and modified RLC elements, respectively. Since there Is no gradient 
in the stresses In the x-dlrectlon, the accuracy Is Independent of the refine- 
ment In the x-dlrectlon. 

As was the case for mechanical load, the assumption of constant Cy 
within an element In the y-dlrectlon severely stiffens the system when only 
one element Is used through the thickness. This Is because the initial gy 
in the upper two lamina is of opposite sign to the Initial gy In the lower 
two lamina. Imposing constant final gy in the y-dlrection results in a 
calculated value of gy = 0, even though physically there is nothing in the 
beam to restrict the expansion or contraction in the y-dlrectlon. Figure 6 
shows that using the modified RLC element eliminates this problem entirely. 
Note that very accurate results are also obtained with just two unmodified RLC 
elements through the thickness. 

The other laminate examined was (H/H/J/J) . This laminate Is like the 
previous one, except the Initial gy is the same in all four plies. For this 
case the assumption of constant final gy in the y-direction Is not a 
problem, since all four laminae are supposed to have nearly the same gy. 
Consequently, even with only an unmodified RLC element through the thickness, 
the error was less than 0.1%. 
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CONCLUSIONS 


A simple two-dimensional (2-D) element was developed for analysis of 
laminated composite materials* The rectangular laminated composite (RLC) 
element eliminates the need to add elements to a model simply to account for 
the material properties of various laminae. Explicit Integration in terms of 
generalized displacements minimizes the algebraic effort required to derive 
the element stliffnesses and the thermal load vector. A substitute shape 
function technique was used to avoid the excessive bending stiffness of 
ordinary bilinear rectangular elements. 

Several sample problems were analyzed using the RLC element. Results 
from these analyses demonstrated that the RLC element accurately accounts for 
the presence of multiple lamina within a single element. Use of the RLC 
element will reduce the number of elements required to analyze many linear 2— D 
laminated composite problems. The basic technique described herein also 
should be applicable for deriving elements for linear three-dimensional (3-D) 
and geometrically nonlinear 2— D and 3-D problems. 
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Table 1 Lamina Properties 


Lamina Type 

Ex 

Jy. 

'’xy 

®xy 

°x 


H 

1.0 

1.0 

.3 

0.385 

l.E-6 

l.E-6 

I 

1.0 

1.0 

.3 

.385 

-l.E-6 

-l.E-6 

J 

1.0 

1.0 

.3 

.385 

-l.E-6 

l.E-6 

S 

.1 

.1 

.3 

.0385 

— 

— 
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(b) 2x 1 mesh 



(c) 4x 1 mesh 



(d) 4 X 2 mesh 



(e) 8 X 1 mesh 


Fig. 4. “ Finite element meshes for .cantilevered 
beam problem. 
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Laminate 
O • (H/H/H/H) 



Mesh 


Fig, 5. “ Accuracy of calculated tip displacement 
for tip loaded cantilevered beam. Open 
s3nnbols are for unmodified RLC element. 
Solid symbols are for modified RLC element. 
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1x1 2x1 4x1 4x2 8x1 
Mesh 


Fig. 6. - Accuracy of calculated tip displacement for 
thermal loading of unsymmetric cantilevered 
beam. Open symbols are for unmodified RLC 
element. Solid symbols are for modified RLC 
element. (Laminate = (H/H/I/I)) 
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